##########################
# quantile estimates
##########################


load("replication_MS/cbq_q1.RData")
q1 = as.data.frame(ms_cbq_stan)
load("replication_MS/cbq_q2.RData")
q2 = as.data.frame(ms_cbq_stan)
load("replication_MS/cbq_q3.RData")
q3 = as.data.frame(ms_cbq_stan)
load("replication_MS/cbq_q4.RData")
q4 = as.data.frame(ms_cbq_stan)
load("replication_MS/cbq_q5.RData")
q5 = as.data.frame(ms_cbq_stan)
load("replication_MS/cbq_q6.RData")
q6 = as.data.frame(ms_cbq_stan)
load("replication_MS/cbq_q7.RData")
q7 = as.data.frame(ms_cbq_stan)
load("replication_MS/cbq_q8.RData")
q8 = as.data.frame(ms_cbq_stan)
load("replication_MS/cbq_q9.RData")
q9 = as.data.frame(ms_cbq_stan)

q1_mean = apply(q1,2,mean)
q2_mean = apply(q2,2,mean)
q3_mean = apply(q3,2,mean)
q4_mean = apply(q4,2,mean)
q5_mean = apply(q5,2,mean)
q6_mean = apply(q6,2,mean)
q7_mean = apply(q7,2,mean)
q8_mean = apply(q8,2,mean)
q9_mean = apply(q9,2,mean)
q_mean = cbind(q1_mean,q2_mean,q3_mean,q4_mean,q5_mean,q6_mean,q7_mean,q8_mean,q9_mean)

q1_25 = apply(q1,2,quantile,0.025)
q2_25 = apply(q2,2,quantile,0.025)
q3_25 = apply(q3,2,quantile,0.025)
q4_25 = apply(q4,2,quantile,0.025)
q5_25 = apply(q5,2,quantile,0.025)
q6_25 = apply(q6,2,quantile,0.025)
q7_25 = apply(q7,2,quantile,0.025)
q8_25 = apply(q8,2,quantile,0.025)
q9_25 = apply(q9,2,quantile,0.025)
q_25 = cbind(q1_25,q2_25,q3_25,q4_25,q5_25,q6_25,q7_25,q8_25,q9_25)

q1_975 = apply(q1,2,quantile,0.975)
q2_975 = apply(q2,2,quantile,0.975)
q3_975 = apply(q3,2,quantile,0.975)
q4_975 = apply(q4,2,quantile,0.975)
q5_975 = apply(q5,2,quantile,0.975)
q6_975 = apply(q6,2,quantile,0.975)
q7_975 = apply(q7,2,quantile,0.975)
q8_975 = apply(q8,2,quantile,0.975)
q9_975 = apply(q9,2,quantile,0.975)
q_975 = cbind(q1_975,q2_975,q3_975,q4_975,q5_975,q6_975,q7_975,q8_975,q9_975)

var_names = c(
  "Minority Coalition",
  "Minimal Winning Coalition",
  "Number of Parties in the Coalition",
  "Largest Party in the Coalition",
  "Median Party in the Coalition",
  "Ideological Divisions in the Coalition",
  "Ideological Divisions within Majority Opposition",
  "Previous Prime Minister in the Coalition",
  "Incumbent Coalition",
  "Minority Coalition where Investiture Vote Required",
  "Anti-System Presence in the Coalition",
  "Pre-Electoral Pact associated with the Coalition",
  "Anti-Pact associated with the Coalition"
)


pdf("figures/figure6_appendix.pdf",width=11,height = 14)
par(mfrow=c(5,3),
    cex.lab=1.5)
plot(NA,xlim=c(0,10),ylim=c(min(q_25[1,]),max(q_975[1,])),
     axes=F,
     xlab="Quantiles",
     ylab=""
     ,main = var_names[1])
axis(1,at=1:9,labels=paste0("0.",1:9))
axis(2)
lines(x=1:9,y=q_mean[1,])
lines(x=1:9,y=q_25[1,],lty="dashed")
lines(x=1:9,y=q_975[1,],lty="dashed")
abline(h=0,lty="dotted",lwd=0.5,col="gray")

plot(NA,xlim=c(0,10),ylim=c(min(q_25[2,]),max(q_975[2,])),
     axes=F,
     xlab="Quantiles",
     ylab=""
     ,main = var_names[2])
axis(1,at=1:9,labels=paste0("0.",1:9))
axis(2)
lines(x=1:9,y=q_mean[2,])
lines(x=1:9,y=q_25[2,],lty="dashed")
lines(x=1:9,y=q_975[2,],lty="dashed")
abline(h=0,lty="dotted",lwd=0.5,col="gray")

plot(NA,xlim=c(0,10),ylim=c(min(q_25[3,]),max(q_975[3,])),
     axes=F,
     xlab="Quantiles",
     ylab=""
     ,main = var_names[3])
axis(1,at=1:9,labels=paste0("0.",1:9))
axis(2)
lines(x=1:9,y=q_mean[3,])
lines(x=1:9,y=q_25[3,],lty="dashed")
lines(x=1:9,y=q_975[3,],lty="dashed")
abline(h=0,lty="dotted",lwd=0.5,col="gray")

plot(NA,xlim=c(0,10),ylim=c(min(q_25[4,]),max(q_975[4,])),
     axes=F,
     xlab="Quantiles",
     ylab=""
     ,main = var_names[4])
axis(1,at=1:9,labels=paste0("0.",1:9))
axis(2)
lines(x=1:9,y=q_mean[4,])
lines(x=1:9,y=q_25[4,],lty="dashed")
lines(x=1:9,y=q_975[4,],lty="dashed")
abline(h=0,lty="dotted",lwd=0.5,col="gray")

plot(NA,xlim=c(0,10),ylim=c(min(q_25[5,]),max(q_975[5,])),
     axes=F,
     xlab="Quantiles",
     ylab=""
     ,main = var_names[5])
axis(1,at=1:9,labels=paste0("0.",1:9))
axis(2)
lines(x=1:9,y=q_mean[5,])
lines(x=1:9,y=q_25[5,],lty="dashed")
lines(x=1:9,y=q_975[5,],lty="dashed")
abline(h=0,lty="dotted",lwd=0.5,col="gray")

plot(NA,xlim=c(0,10),ylim=c(min(q_25[6,]),max(q_975[6,])),
     axes=F,
     xlab="Quantiles",
     ylab=""
     ,main = var_names[6])
axis(1,at=1:9,labels=paste0("0.",1:9))
axis(2)
lines(x=1:9,y=q_mean[6,])
lines(x=1:9,y=q_25[6,],lty="dashed")
lines(x=1:9,y=q_975[6,],lty="dashed")
abline(h=0,lty="dotted",lwd=0.5,col="gray")

plot(NA,xlim=c(0,10),ylim=c(min(q_25[7,]),max(q_975[7,])),
     axes=F,
     xlab="Quantiles",
     ylab=""
     ,main = var_names[7])
axis(1,at=1:9,labels=paste0("0.",1:9))
axis(2)
lines(x=1:9,y=q_mean[7,])
lines(x=1:9,y=q_25[7,],lty="dashed")
lines(x=1:9,y=q_975[7,],lty="dashed")
abline(h=0,lty="dotted",lwd=0.5,col="gray")

plot(NA,xlim=c(0,10),ylim=c(min(q_25[8,]),max(q_975[8,])),
     axes=F,
     xlab="Quantiles",
     ylab=""
     ,main = var_names[8])
axis(1,at=1:9,labels=paste0("0.",1:9))
axis(2)
lines(x=1:9,y=q_mean[8,])
lines(x=1:9,y=q_25[8,],lty="dashed")
lines(x=1:9,y=q_975[8,],lty="dashed")
abline(h=0,lty="dotted",lwd=0.5,col="gray")

plot(NA,xlim=c(0,10),ylim=c(min(q_25[9,]),max(q_975[9,])),
     axes=F,
     xlab="Quantiles",
     ylab=""
     ,main = var_names[9])
axis(1,at=1:9,labels=paste0("0.",1:9))
axis(2)
lines(x=1:9,y=q_mean[9,])
lines(x=1:9,y=q_25[9,],lty="dashed")
lines(x=1:9,y=q_975[9,],lty="dashed")
abline(h=0,lty="dotted",lwd=0.5,col="gray")

plot(NA,xlim=c(0,10),ylim=c(min(q_25[10,]),max(q_975[10,])),
     axes=F,
     xlab="Quantiles",
     ylab=""
     ,main = var_names[10])
axis(1,at=1:9,labels=paste0("0.",1:9))
axis(2)
lines(x=1:9,y=q_mean[10,])
lines(x=1:9,y=q_25[10,],lty="dashed")
lines(x=1:9,y=q_975[10,],lty="dashed")
abline(h=0,lty="dotted",lwd=0.5,col="gray")

plot(NA,xlim=c(0,10),ylim=c(min(q_25[11,]),max(q_975[11,])),
     axes=F,
     xlab="Quantiles",
     ylab=""
     ,main = var_names[11])
axis(1,at=1:9,labels=paste0("0.",1:9))
axis(2)
lines(x=1:9,y=q_mean[11,])
lines(x=1:9,y=q_25[11,],lty="dashed")
lines(x=1:9,y=q_975[11,],lty="dashed")
abline(h=0,lty="dotted",lwd=0.5,col="gray")

plot(NA,xlim=c(0,10),ylim=c(min(q_25[12,]),max(q_975[12,])),
     axes=F,
     xlab="Quantiles",
     ylab=""
     ,main = var_names[12])
axis(1,at=1:9,labels=paste0("0.",1:9))
axis(2)
lines(x=1:9,y=q_mean[12,])
lines(x=1:9,y=q_25[12,],lty="dashed")
lines(x=1:9,y=q_975[12,],lty="dashed")
abline(h=0,lty="dotted",lwd=0.5,col="gray")

plot(NA,xlim=c(0,10),ylim=c(min(q_25[13,]),max(q_975[13,])),
     axes=F,
     xlab="Quantiles",
     ylab=""
     ,main = var_names[13])
axis(1,at=1:9,labels=paste0("0.",1:9))
axis(2)
lines(x=1:9,y=q_mean[13,])
lines(x=1:9,y=q_25[13,],lty="dashed")
lines(x=1:9,y=q_975[13,],lty="dashed")
abline(h=0,lty="dotted",lwd=0.5,col="gray")

plot(NA,xlim=c(0,10),ylim=c(0,10),
     axes=F,
     xlab="",
     ylab="")

legend("topleft",legend=c("Posterior Mean","95% Credible Interval"),
       lty=c("solid","dashed"),
       bty="n",
       cex=1.4)

dev.off()







